home *** CD-ROM | disk | FTP | other *** search
/ Atari Mega Archive 1 / Atari Mega Archive - Volume 1.iso / mint / lib / mntlib44.zoo / mntlib / atof.c < prev    next >
C/C++ Source or Header  |  1993-09-15  |  17KB  |  793 lines
  1. /*
  2.  *     double strtod (str, endptr);
  3.  *     const char *str;
  4.  *     char **endptr;
  5.  *        if !NULL, on return, points to char in str where conv. stopped
  6.  *
  7.  *     double atof (str)
  8.  *     const char *str;
  9.  *
  10.  * recognizes:
  11.  [spaces] [sign] digits [ [.] [ [digits] [ [e|E|d|D] [space|sign] [int][F|L]]]]
  12.  *
  13.  * returns:
  14.  *    the number
  15.  *        on overflow: HUGE_VAL and errno = ERANGE
  16.  *        on underflow: -HUGE_VAL and errno = ERANGE
  17.  *
  18.  * bugs:
  19.  *   naive over/underflow detection
  20.  *
  21.  *    ++jrb bammi@dsrgsun.ces.cwru.edu
  22.  *
  23.  * 07/29/89:
  24.  *    ok, before you beat me over the head, here is a hopefully
  25.  *    much improved one, with proper regard for rounding/precision etc
  26.  *    (had to read up on everything i did'nt want to know in K. V2)
  27.  *    the old naive coding is at the end bracketed by ifdef __OLD__.
  28.  *    modeled after peter housels posting on comp.os.minix.
  29.  *    thanks peter!
  30.  *
  31.  *    mjr: 68881 version added
  32.  *
  33.  *    schwab: 68881 version replaced with new implementation that
  34.  *        uses fmovep for maximum precision, no bits lost anymore!
  35.  */
  36.  
  37. #if !defined (__M68881__) && !defined (sfp004)
  38.  
  39. #if !(defined(unix) || defined(minix))
  40. #include <stddef.h>
  41. #include <stdlib.h>
  42. #include <float.h>
  43. #endif
  44. #include <errno.h>
  45. #include <assert.h>
  46. #include <math.h>
  47.  
  48. #ifdef minix
  49. #include "lib.h"
  50. #endif
  51.  
  52. #ifndef _COMPILER_H
  53. #include <compiler.h>
  54. #endif
  55.  
  56. extern int errno;
  57. #if (defined(unix) || defined(minix))
  58. #define NULL     ((void *)0)
  59. #endif
  60.  
  61. #define Ise(c)        ((c == 'e') || (c == 'E') || (c == 'd') || (c == 'D'))
  62. #define Isdigit(c)    ((c <= '9') && (c >= '0'))
  63. #define Isspace(c)    ((c == ' ') || (c == '\t'))
  64. #define Issign(c)    ((c == '-') || (c == '+'))
  65. #define IsValidTrail(c) ((c == 'F') || (c == 'L'))
  66. #define Val(c)        ((c - '0'))
  67.  
  68. #define MAXDOUBLE    DBL_MAX
  69. #define MINDOUBLE    DBL_MIN
  70.  
  71. #define MAXF  1.797693134862316
  72. #define MINF  2.225073858507201
  73. #define MAXE  308
  74. #define MINE  (-308)
  75.  
  76. /* another alias for ieee double */
  77. struct ldouble {
  78.     unsigned long hi, lo;
  79. };
  80.  
  81. static int __ten_mul __PROTO((double *acc, int digit));
  82. static double __adjust __PROTO((double *acc, int dexp, int sign));
  83.  
  84. #ifdef __OLD__
  85. static double __ten_pow __PROTO((double r, int e));
  86. #endif
  87.  
  88. /*
  89.  * mul 64 bit accumulator by 10 and add digit
  90.  * algorithm:
  91.  *    10x = 2( 4x + x ) == ( x<<2 + x) << 1
  92.  *    result = 10x + digit
  93.  */
  94. static int __ten_mul(acc, digit)
  95. double *acc;
  96. int digit;
  97. {
  98.     register unsigned long d0, d1, d2, d3;
  99.     register          short i;
  100.     
  101.     d2 = d0 = ((struct ldouble *)acc)->hi;
  102.     d3 = d1 = ((struct ldouble *)acc)->lo;
  103.  
  104.     /* check possibility of overflow */
  105. /*    if( (d0 & 0x0FFF0000L) >= 0x0ccc0000L ) */
  106. /*    if( (d0 & 0x70000000L) != 0 ) */
  107.     if( (d0 & 0xF0000000L) != 0 )
  108.     /* report overflow somehow */
  109.     return 1;
  110.     
  111.     /* 10acc == 2(4acc + acc) */
  112.     for(i = 0; i < 2; i++)
  113.     {  /* 4acc == ((acc) << 2) */
  114.     asm volatile("    lsll    #1,%1;
  115.              roxll    #1,%0"    /* shift L 64 bit acc 1bit */
  116.         : "=d" (d0), "=d" (d1)
  117.         : "0"  (d0), "1"  (d1) );
  118.     }
  119.  
  120.     /* 4acc + acc */
  121.     asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  122.     asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  123.  
  124.     /* (4acc + acc) << 1 */
  125.     asm volatile("  lsll    #1,%1;
  126.              roxll   #1,%0"    /* shift L 64 bit acc 1bit */
  127.         : "=d" (d0), "=d" (d1)
  128.         : "0"  (d0), "1"  (d1) );
  129.  
  130.     /* add in digit */
  131.     d2 = 0;
  132.     d3 = digit;
  133.     asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  134.     asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  135.  
  136.  
  137.     /* stuff result back into acc */
  138.     ((struct ldouble *)acc)->hi = d0;
  139.     ((struct ldouble *)acc)->lo = d1;
  140.  
  141.     return 0;    /* no overflow */
  142. }
  143.  
  144. #include "flonum.h"
  145.  
  146. static double __adjust(acc, dexp, sign)
  147. double *acc;    /* the 64 bit accumulator */
  148. int     dexp;    /* decimal exponent       */
  149. int    sign;    /* sign flag          */
  150. {
  151.     register unsigned long  d0, d1, d2, d3;
  152.     register          short i;
  153.     register           short bexp = 0; /* binary exponent */
  154.     unsigned short tmp[4];
  155.     double r;
  156.  
  157. #ifdef __STDC__
  158.     double __normdf( double d, int exp, int sign, int rbits );
  159.     double ldexp(double d, int exp);
  160. #else
  161.     extern double __normdf();
  162.     extern double ldexp();
  163. #endif    
  164.     d0 = ((struct ldouble *)acc)->hi;
  165.     d1 = ((struct ldouble *)acc)->lo;
  166.     while(dexp != 0)
  167.     {    /* something to do */
  168.     if(dexp > 0)
  169.     { /* need to scale up by mul */
  170.         while(d0 > 429496729 ) /* 2**31 / 5 */
  171.         {    /* possibility of overflow, div/2 */
  172.         asm volatile(" lsrl    #1,%1;
  173.                     roxrl    #1,%0"
  174.             : "=d" (d1), "=d" (d0)
  175.             : "0"  (d1), "1"  (d0));
  176.         bexp++;
  177.         }
  178.         /* acc * 10 == 2(4acc + acc) */
  179.         d2 = d0;
  180.         d3 = d1;
  181.         for(i = 0; i < 2; i++)
  182.         {  /* 4acc == ((acc) << 2) */
  183.         asm volatile("    lsll    #1,%1;
  184.                  roxll    #1,%0"    /* shift L 64 bit acc 1bit */
  185.                  : "=d" (d0), "=d" (d1)
  186.                  : "0"  (d0), "1"  (d1) );
  187.         }
  188.  
  189.         /* 4acc + acc */
  190.         asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  191.         asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  192.  
  193.         /* (4acc + acc) << 1 */
  194.         bexp++;    /* increment exponent to effectively acc * 10 */
  195.         dexp--;
  196.     }
  197.     else /* (dexp < 0) */
  198.     { /* scale down by 10 */
  199.         while((d0 & 0xC0000000L) == 0)
  200.         {    /* preserve prec by ensuring upper bits are set before div */
  201.         asm volatile(" lsll  #1,%1;
  202.                     roxll #1,%0" /* shift L to move bits up */
  203.             : "=d" (d0), "=d" (d1)
  204.             : "0"  (d0), "1"  (d1) );
  205.         bexp--;    /* compensate for the shift */
  206.         }
  207.         /* acc/10 = acc/5/2 */
  208.         *((unsigned long *)&tmp[0]) = d0;
  209.         *((unsigned long *)&tmp[2]) = d1;
  210.         d2 = (unsigned long)tmp[0];
  211.         asm volatile(" divu #5,%0" : "=d" (d2) : "0" (d2));
  212.         tmp[0] = (unsigned short)d2;    /* the quotient only */
  213.         for(i = 1; i < 4; i++)
  214.         {
  215.         d2 = (d2 & 0xFFFF0000L) | (unsigned long)tmp[i]; /* REM|next */
  216.         asm volatile(" divu #5,%0" : "=d" (d2) : "0" (d2));
  217.         tmp[i] = (unsigned short)d2;
  218.         }
  219.         d0 = *((unsigned long *)&tmp[0]);
  220.         d1 = *((unsigned long *)&tmp[2]);
  221.         /* acc/2 */
  222.         bexp--;    /* div/2 taken care of by decrementing binary exp */
  223.         dexp++;
  224.     }
  225.     }
  226.     
  227.     /* stuff the result back into acc */
  228.     ((struct ldouble *)acc)->hi = d0;
  229.     ((struct ldouble *)acc)->lo = d1;
  230.  
  231.     /* normalize it */
  232.     r = __normdf( *acc, ((0x03ff - 1) + 53), (sign)? -1 : 0, 0 );
  233.     /* now shove in the binary exponent */
  234.     return ldexp(r, bexp);
  235. }
  236.  
  237. /* flags */
  238. #define SIGN    0x01
  239. #define ESIGN    0x02
  240. #define DECP    0x04
  241. #define CONVF    0x08
  242.  
  243. double strtod (s, endptr)
  244. register const char *s;
  245. char **endptr;
  246. {
  247.     double         accum = 0.0;
  248.     register short flags = 0;
  249.     register short texp  = 0;
  250.     register short e     = 0;
  251.     double       zero = 0.0;
  252.     const char        *tmp;
  253.  
  254.     assert ((s != NULL));
  255.  
  256.     if(endptr != NULL) *endptr = (char *)s;
  257.     while(Isspace(*s)) s++;
  258.     if(*s == '\0')
  259.     {    /* just leading spaces */
  260.     return zero;
  261.     }
  262.  
  263.     if(Issign(*s))
  264.     {
  265.     if(*s == '-') flags = SIGN;
  266.     if(*++s == '\0')
  267.     {   /* "+|-" : should be an error ? */
  268.         return zero;
  269.     }
  270.     }
  271.  
  272.     do {
  273.     if (Isdigit(*s))
  274.     {
  275.         flags |= CONVF;
  276.         if( __ten_mul(&accum, Val(*s)) ) texp++;
  277.         if(flags & DECP) texp--;
  278.     }
  279.     else if(*s == '.')
  280.     {
  281.         if (flags & DECP)  /* second decimal point */
  282.         break;
  283.         flags |= DECP;
  284.     }
  285.     else
  286.         break;
  287.     s++;
  288.     } while (1);
  289.  
  290.     if(Ise(*s))
  291.     {
  292.     if(*++s != '\0') /* skip e|E|d|D */
  293.     {  /* ! ([s]xxx[.[yyy]]e)  */
  294.          tmp = s;
  295.         while(Isspace(*s)) s++; /* Ansi allows spaces after e */
  296.         if(*s != '\0')
  297.         { /*  ! ([s]xxx[.[yyy]]e[space])  */
  298.  
  299.         if(Issign(*s))
  300.             if(*s++ == '-') flags |= ESIGN;
  301.  
  302.         if(*s != '\0')
  303.         { /*  ! ([s]xxx[.[yyy]]e[s])  -- error?? */
  304.  
  305.             for(; Isdigit(*s); s++)
  306.             e = (((e << 2) + e) << 1) + Val(*s);
  307.  
  308.             if(IsValidTrail(*s)) s++;
  309.         
  310.             /* dont care what comes after this */
  311.             if(flags & ESIGN)
  312.             texp -= e;
  313.             else
  314.             texp += e;
  315.         }
  316.         }
  317.         if(s == tmp) s++;    /* back up pointer for a trailing e|E|d|D */
  318.     }
  319.     }
  320.  
  321.     if((endptr != NULL) && (flags & CONVF))
  322.     *endptr = (char *) s;
  323.     if(accum == zero) return zero;
  324.  
  325.     return __adjust(&accum, (int)texp, (int)(flags & SIGN));
  326. }
  327.  
  328. double atof(s)
  329. const char *s;
  330. {
  331. #ifdef __OLD__
  332.     extern double strtod();
  333. #endif
  334.     return strt